Master Equation Approach to Molecular Motors 
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A master equation approach to molecular motors allows to describe a mechano-chemical cyclic 
system where chemical and translational degrees of freedom are treated on an equal footing. A 
generalized detailed balance condition in the out of equilibrium regime is shown to be compatible 
with the Fokker-Planck equation in the continuum limit. The Onsager reciprocity relations hold 
for stationary states close to equilibrium, provided the generalized detailed balance condition is 
satisfied. Semi-phenomenological considerations in the case of motor proteins lead to a discrete 
kinetics model, for which interesting observable quantities can be directly calculated and compared 
with experimental data. 



PACS numbers: 87.16.Nn, 87.15.Aa, 05.10.Gg, 05 

Current models used to describe the properties of 
molecular motors and the energy transduction process 
fall under two distinct categories: continuous mod- 
els |l]-|3| and discrete models pLEj. Both represent a 
coarse grained description of a very complicated physico- 
chemical system and the use of one or the other depends 
on the quantities one is interested in. For example, con- 
tinuous models are very useful to investigate the role of 
an external force on chemical kinetics Q , since the exter- 
nal force is inserted into the Fokker-Planck (FP) equa- 
tions without any ambiguity. This is no longer true for 
discrete models, when one has to resort to some ad hoc 
principle or a priori reasoning to insert force in transi- 
tion rates p|. Nonetheless discrete models present the 
important advantage of being analytically solvable, as it 
happens, for instance, in jump processes Q. On the other 
hand an analytical solution is quite difficult to obtain in 
the general case of continuous models, and one has to 
resort to complex numerical integrations. In this paper 
we introduce a discrete model, similar to the ones pro- 
posed in Q and ||, but with the following constraint: if 
a is the lattice distance between subsequent spatial po- 
sitions of the system, the continuous model should be 
obtained as a limit of the discrete one for a — > 0. The 
connection between a kinetic theory involving activated 
transitions over potential energy barriers and a diffusion 
theory approach based on a FP equation dates back to 
Kramers see || for a recent review. In M], the idea 
was applied to models for protein motors, but the force 
dependence was left in equal apportionments over back- 
ward and forward transition rates and the experimental 
results on the force dependence of the apparent Michaelis 
constant [[ll| were not available. In JL2|| , a general the- 
ory for motor proteins was presented. This theory was 
developed in a two dimensional manifold and complex 
integrations over state variables were used to calculate 
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force dependent transition rates over potential barriers 
in a discrete model. In the present paper we do not use 
complex integrations; instead, we identify in the general- 
ized detailed balance the condition for a discrete model 
to be compatible with a continuous one. Few parameters 
arc needed to capture the overall shape of the potential 
energy surface which is no longer needed in full detail. 
Detailed balance was used in jl3| to calculate transition 
rates for particles diffusing over potential barriers. The 
obtained discrete kinetics model led to a fast and reli- 
able numerical procedure to calculate mean velocity of 
correlation ratchets. A model similar to ours was also 
developed in [|l4| in the context of thermal ratchets, even 
if no connection with the actual chemistry of motor pro- 
teins was done. Moreover, motor proteins are isother- 
mal, and therefore arc better described by correlation 
ratchets. In |l5| , a master equation approach was used 
to investigate the force generation in RNA polymerase, 
which can be considered a motor protein, even if it differs 
from kinesin, myosin and dynein both in structure and 
function. In section | we outline the general framework, 
which may be useful not only for modeling molecular 
motors, but also any mechano-chemical cyclic system. 
In the general formulation of our model, the chemical re- 
action coordinate is treated on an equal footing as the 
spatial one. Onsager reciprocity relations will be shown 
to hold in the most general case in the stationary peri- 
odic close to equilibrium state, provided detailed balance 
is verified. In section |n| we specify our model to the con- 
text of molecular motors. We show that our model may 
be regarded as the discrete analogue of the continuous 
one proposed in ^|J^], leading to a clear interpretation 
of the generalized forces and currents introduced in sec- 
tion H In section III a discrete chemical kinetics model 
with a generalized detailed balance condition is defined 
for the case of motor proteins. Semi-phenomenological 
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considerations help to decide the apportionments of gen- 
eralized forces over forward and backward transitions. In 
section |y] two example models are studied and their pre- 
dictions compared with experimental data. 



I. GENERAL FRAMEWORK 

Our model will describe the time evolution of a thermo- 
dynamic out-of-equilibrium system in a complex phase 
space. The state of the system (a motor protein, an ion 
pump or whatsoever) is determined by the thermody- 
namic parameters x a (a = l,...,d). One of these pa- 
rameters may represent the position of the protein cen- 
ter of mass, another a chemical reaction coordinate (or 
a variable indicating the conformation of the protein) 
and so on. In general the number of these thermody- 
namic parameters (and hence the dimensionality of the 
system under study) is sufficient to identify the state of 
the system by a direct experimental measure. Therefore 
we will assume the state of the system to be described 
by a d— dimensional vector X subject to a time evolu- 
tion in a dimensional discrete phase space, which can 
be mapped on Z d (in general the lattice spacing in each 
direction will be different). 

The probability of being in a particular state-X at 
time t is written as Px(t). Since the system must be in 
one of the X states, the normalization condition follows: 



E^) = 1 



V/ 



(1) 



Wxy is defined as the transition probability per unit 
time from state Y to state X and it is assumed to be 
time independent. 

Another hypothesis is the full periodicity along any di- 
rection a. This is usually assumed for all models of motor 
proteins, (see jl0],[l2]] ) . The periodicity N a depends on a, 
but we assume N a > 1, since it is always possible to re- 
duce the steps until this constraint is satisfied. In other 
words we assume that the state described by the parame- 
ters (/ii\Ti+a;i,...,/ <i iV d + x <i ) with L = {h,...l d } G Z d 
is equivalent to the state described by the parameters 
(xi, . . .,Xd). 

We introduce the variable: 



1 for 1 < x a < N a Vq G {1, . . . , d} 
otherwise 



(2) 



which is an indicator of the period in which the system 
is moving and will be useful for subsequent calculations. 

The time evolution of the system is simply given by 
the master equation: 

Px = J2 Y ( w xyPy - WyxPx) = E Y l xyPy, (3) 
where we have defined: 



L X y = Wxy ~ S X y E z W Z x- 
From this definition, it follows: 

E^ = °^E A -^ = °> 



(4) 



(5) 



which is consistent with the normalization condition, 
cq. ([!]). Since the system is assumed to be periodic, 



Wx+ln,y+ln = Wxy VX g Z° 



(6) 



where LN = (hNi, . . . , IdNd). It is simple to show that 
the same rule holds also for the matrix Lxy, defined in 
eq. (§. 

We introduce a time independent variable qx , depend- 
ing explicitly on the X coordinate, i.e. on the state of 
the system. We define the current J q conjugated to the 
variable qx- 



Jq ~ dt ' 



(7) 



where (q) is the average: (q) = ^2 x q X Px- By applying 
eqs. (||) and (0) it easily follows: 



XY^Y ■ 



(8) 



We introduce probabilities and transition rates over all 
periods, following some of the formalism of the d — 1 case 
studied in : 



fix, VP 



X+LN 



£-XY — E L Lx,Y+LN- 



(9) 
(10) 



By definition, Rx is a periodic quantity. It is easy 
to show that also Cxy is periodic in both arguments X 
and Y . At variance of Lxy, C-xy is a finite matrix; also 
Rx is a finite vector, whereas Px is not. From the time 
evolution of Px, we can easily obtain the time evolution 
of R x : 



Rx — Ey LxvRi 



For any variable fy and using eq. 
property holds: 



E x $ x E x Xx E L f- 



X+LN- 



(11) 

the following 
(12) 



Applying property (|T^) to eq. (^Tj) and the periodicity 
of R, we obtain: 



Rx — Ev ^ XY ^ Y ' 



(13) 



where, by definition, ^2y ~ X Y ■ i- c - a primed sum 
is restricted only to one period along any axis. This 
is a master equation for a system with a finite number 
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(n 



a=l 



N a ) of states. The Cxy matrix is finite and has 



the following properties: 



E 

^x 



C XY > iorX ^ Y 
Cxy = VY~ e Z d . 



It is easy to show, by applying eqs. (12) and 



E' ^ = E 



p 



X 



1. 



(14) 
(15) 

that: 
(16) 



Therefore, there exists a stationary solution Rx of 
cq.(^) and, under general hypotheses (always satis- 
fied in the examples treated in the next sections), it is 
unique Jl(| . It is also periodic, by definition. 

If qx = x a for any value of a, then qx~QY = Qx+ln — 
Qy+ln is also periodic, and we can rewrite the current 
as: 

J q = ^Z XY x Y (Qx -qy)LxYRY , qx = x a . (17) 

A subsequent application of property ( |l2| ) to eq. (ft 

gives: 



Jg = E xy Xx ( qx ~ qY ^ LxyRy 



qx 



(18) 



and, after some manipulation, we obtain the following 
expression: 



Jo 



Xx + Xy 



^XY 



{qx - qy) {LxyRy 



LyxRx) , 
(19) 



where the argument in the sum is evidently symmetric 
under a change X Y. 

If the detailed balance condition for the periodic sta- 
tionary state, Rx = Rx, holds: 



WxyRy = W yx Rx, 
which is equivalent to the following: 

LxyRy = LyxRx, 



(20) 



(21) 



then from eq. (|T^) the net stationary current is zercj]. 
A net flow, i.e. a nonzero stationary current, may oc- 
cur only if the detailed balance condition, eq. (p0|), is 
violated. This can be done in several ways: in the con- 
tinuous model proposed in rcf. ||, for instance, detailed 
balance holds separately for each chemical reaction, in- 
troducing the chemical potential AyU. In our model we 
introduce a set of generalized forces, able to drive the sys- 
tem out of equilibrium, so that a finite stationary current 



may occur. Each generalized force, f a is coupled to one 
generalized coordinate x a . Both the transition matrix 
Lxy and the stationary solution will depend explicitly 
on the force vector F, so that: 



£xy(F)R y (F) = 



(22) 



Of course, at equilibrium F — and the stationary 
currents are all identically zero. Our assumption is that 
condition is replaced by a generalized detailed bal- 
ance condition (in this section and in appendix the 
factor (3 = I/UbT is absorbed in the definition of F): 



L XY (F)R Y (0y y = L YX (F)R x (0)e 



FX 



(23) 



We remark that the stationary solution Rx{F) does 
not satisfy a detailed balance condition: 



L xy {F)Ry{F) + L YX {F)Rx{F) 



(24) 



An equality in eq. (gj) would imply that Ry{F) could 
be written as: 



Rx{F) cxi? x (0)e 



FX 



(25) 



which is, evidently, not periodic. 

In a local thermodynamic equilibrium the generalized 
stationary currents can be written in the following form: 



Jx a — E ^apfp- 

13=1 



(26) 



The coefficients \ a p are called Onsagcr coefficients. In 
general, even when the linear approximation cannot be 
applied, we can define: 



A 



dJ x 



a/3 



(27) 



We show in appendix^ that, provided the generalized 
detailed balance condition, eq. (^3|), holds, these coeffi- 
cients verify the generalized Gyarmati-Li JT^ ] reciprocal 
relations: 



Ktp = A^c Va, f3 G {1, . . . ,d} 



(28) 



and, hence, the Onsager reciprocity relations. These 
properties are general and do not depend on the specific 
parameters of the model and are based on the general- 
ized detailed balance condition eq. ((2^). This condition 
is a common assumption also for continuous models, as 
discussed further in the next section. 



*Notice that if eq. ( pl[ ) holds, then also CxyRy = LyxRx 
holds, but the converse is not guaranteed to be true. 
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II. CONTINUUM LIMIT FOR MOLECULAR 
MOTORS 

The usual choice for molecular motors is a 2-dimen- 
sional manifold in which one direction represents the po- 
sition of the center of mass along the linear track (the 
microtubule or the actin filament). The other is the reac- 
tion coordinate for the ATP hydrolysis, (see fll2)), which 
is also related to the conformational changes of the mo- 
tor protein. These conformational changes are commonly 
thought to occur after binding ATP and release of reac- 
tion products ADP and Pi @. 

The transition rates will be therefore written as W*"™ 
where the subscript xy denotes a transition from spatial 
position y to spatial position x whereas the superscript 
nm stands for a transition from a state with m ATP 
molecules to a state with n ones^]. Of course this vari- 
able may also represent a non-integer chemical reaction 
coordinate (accounting for multiple states models), but 
this is the simplest possible choice. The periodicity is 
not specified for the spatial direction, while it is 1 for the 
chemical direction, i.e. we are assuming that the state of 
the motor in presence of n ATP molecules is equivalent to 
the one in presence of n + 1 ATP molecules. We remark 
that this assumption does not mean that a thermody- 
namic system with n ATP molecules is equivalent to the 
same thermodynamic system with n + 1 ATP molecules, 
but only that the chemical state of the motor protein af- 
ter the reaction cycle is completed (and 1 ATP molecule 
is consumed or produced) is equivalent to the state it 
was before entering the cycle. We allow only transitions 
from a position x to x + a and x — a. All other transi- 
tion rates will be identically 0. Chemical transitions, i.e. 
those involving the superscripts nm (with n =/= m) will 
be specified later. At the moment we only use transition 
rates in the forrrffl w xy = Y^m ■ According to the 
definition eq. (^|) and due to our choice of periodicity 1 
in the chemical reaction coordinate, the master equation 
for the rate of change of R x , or cq. (JTT|) is: 



(29) 



R x — 'Wx^x+aRx+a ^x-\-a,xRx ~f~ 

~t~ ^x.x — aRx — a ^x~-a,xRx 

We define a discrete current j x : 

jx {lVx-\-a,xRx ^a;,a;+a-Ra;+a) j 



(30) 



Applying this definition to cq. (E9| 



(jx ~ jx-a) = -a(Vdj)x 



(31) 



where (Vdj)x is by definition the discrete gradient of j x . 
In the continuum (FP) description: 

P(x,t) = -Vj(x,t), (32) 
j(x, t) = ^ [-TWP(x, t) - P(x, t)VV{x) + P{x, t)f] 



T 



(33) 



where T is the temperature, D the diffusion constant, 
V(x) a periodic potential, / the force and j(x,t) is the 
probability density current. The potential V(x) has the 
same period as the transition rates w xy . 
For very small values of a we define: 



R x = aP(x, t) 



(34) 



If in cq. (|30|) we expand R x + a to the first order in a, 
the discrete current can be rewritten as: 

jx = Rx(wx+ a ,x - w x ^ x+a ) - a(VR x )w x , x+a (35) 



Using eqs. (|3l]), (|3^) and (p4[), for very small a, one 
gets: 



■ j(x,t). 



(36) 



By means of this correspondence, eq. (pq), and using 
eqs. (p3|) and (pa), we obtain: 



D = a w x . x+a 



w 



x.x-\- 



a) 



-DVV(x) +Df = a (w x+a 
Solving for w x+a , x : 

w x + a ,x = w x . x+a [l - a/3(VV(x) - /)], 



(37) 
(38) 



(39) 



where eq. ( p7\ ) has been used. If we suppose that the 
following relation holds: 



-0(V(x+a)-V(x)-af) 



(40) 



then, eq. ( |39| ) is satisfied in the continuum limit (a — > 0). 

This relation corresponds to the standard detailed bal- 
ance condition when / = and the stationary periodic 
solution is: 



Rx{0) 



-/3V{x 



(41) 



whereas eq. ( fl(]| ) corresponds to a generalized detailed 
balance condition, eq. (|23|), when / ^ 0. 



T Using the definitions of the previous section the spatial 
and chemical directions correspond to a — 1 and ot = 2 
respectively. 

"''The dependence on n is no more necessary, since W£™ is 
periodic in n and m with period 1. 
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A. The chemical reaction 



By analogy with the continuum model M , we suppose 
that each transition rate W£™ is essentially due to 3 sub- 
processes, leading to ATP consumption (a transitions in 
||), ATP production (7) and no change in ATP concen- 
tration or thermal transitions (/?). 

Since a chemical reaction is present in a and 7 pro- 
cesses, the chemical potential difference, Afi — (J-atp — 
{J-adp — A* Pj, plays the role of a generalized force, con- 
jugated to the number of ATP molecules. In this paper 
only transitions from m = n + An to n ATP molecules 
with An = ±1,0 will be considered. 

The generalized detailed balance condition eq. ( |23| ) is 
therefore: 



x-\-a,x 



w. 



n+An,n 
x,x-\-a 



-/3(V(x+a)-V(x)-fa~AnAfi) 



(42) 



where e~^ v ^ oc R x (0) is the stationary equilibrium 
solution when all generalized forces are zero. Transi- 
tions with An > 0(< 0) correspond to ATP consump- 
tion (production). Eq. ( ^2[ ) states that when An > 
(A/z < 0) transitions leading to ATP consumption (pro- 
duction) are more favorable and lead to a spatial ad- 
vancement of the motor. This is the core of the energy 
transduction process: chemical energy is used to perform 
mechanical work against a load —/or chemical energy is 
produced performing a mechanical work on the protein. 
Notice that, according to our choice of notation, fi = (3f 
and /2 = — f3 An. 

These generalized detailed balance conditions are the 
same introduced in j3j, in the limit a — > 0. We remark 
that this scheme corresponds to a periodicity 1 in the 
chemical reaction coordinate. A different periodicity in 
the chemical coordinate would allow to introduce differ- 
ent potential shapes V(x, £) for each value of the reaction 
coordinate £, as for instance in [^). Using periodicity 2 in 
the chemical coordinate, the continuous two-state model 
[^|,D is fully recovered in the a — > limit with the same 
detailed balance conditions as in PL The general 



^Indeed, using periodicity 2 for the chemical coordinate, if 
we allow only transitions between neighboring sites both in 
the spatial and chemical coordinate, eq. (jl^) becomes: 

R 1 = W 11 . R 1 ^ - W 1 ! R 1 +W 11 R 1 + 

v x * ' x,X-\-a J - b a?+a r ' x + a^x^^x 1 rr x,x — a^^x — a 1 

- Wl\ x Rl + WllRl - WllR\ (43) 



E>2 __ T,T/22 t-,2 TT/22 p2 . tt^22 t->2 , 

1 ' ' x ,x-\-a- r *'x-\-a ' ' x-^a^x-^x > vv x,x — a^x — a ' 

- W™Rl + WllRl - WllRl 



(44) 



where 1 and 2 represent the two conformational state indexes 
and not the number of ATP molecules. In the continuum 
limit, this model is perfectly equivalent to the ones proposed 



scheme introduced here allows to treat the chemical and 
mechanical coordinates on an equal footing and write the 
detailed balance condition in a more standard and trans- 
parent way. 



B. Generalized currents 

According to the definitions given in section |, the sta- 
tionary current J t (t stays for translational motion), as- 
sociated to the protein center of mass corresponds to the 
velocity of the motor protein. Its full expression is given 
by: 



Jt = a 



[UJx-\-a,x 



R, 



(45) 



whereas the stationary current J c associated to the chem- 
ical coordinate is: 



Jr. = 



W n,n+1 w n,n+l 
r r x-\-a,x 1 r r x — a,x 



w. 



x-\-a,x 



w: 



R., 



(46) 



vv xy 



which does not depend on n since W£™ 

This current corresponds to the number of ATP mole- 
cules consumed per unit time, so we will refer to it as the 
"rate of ATP consumption" . Notice that this current has 
the opposite sign of the one in eq.(19). This is consistent 
with the above choice /2 = — /3A/x. We remark that /3 
transitions do not contribute to ATP production since 
they do not involve any chemical reaction. These consid- 
erations will be used in the next section to develop and 
fully characterize a very simple discrete model whose con- 
tinuum limit is still described by a Fokker-Planck equa- 
tion. 



III. DEFINITION OF THE MODEL 

We showed that eq. (ppj) , i.e. a generalized detailed 
balance condition, is sufficient for a discrete model to be 
compatible with a Fokker Planck equation in the contin- 
uum limit. Actually eq. ( f40| ) is compatible with a very 
general class of models, for which any apportionment of 
force / and chemical potential difference Au over for- 
ward (in space) and backward transitions, is perfectly 
reasonable. This is true as long as each transition al- 
lows the protein to take a substep which can be made 
infinitcsimally small. If this is not possible for some of 



with the following substitutions: w x ,x±a 



V(x) = Vi(x), i = 1 

= Wi(x) in eqs 
The effective mobility £ 



2 in eq. @ and Wl 

in the notation of refs 
of refs. Ill is f5a 2 Wl] x+a . 
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these substeps, the simple Fokker-Planck equation may 
not be able to describe the system in the continuum limit. 
Nevertheless the Onsager relations still hold, as long as 
cq. (^3|) holds and provided that the sum of all substeps 
is equal to one period. The continuum limit, left alone, 
is therefore not sufficient to fix the parameters required 
to define discrete models. 

To be as general as possible, we use forward (uj) and 
backward (wj) transition rates, with a priori not specified 
apportionments: 

Uj (f, A/i) = Wj(f, ArfeW-i^^N (47) 
Wj (f, A M ) - wj(f, A»)e^- a 7f- m 7 A »\ (48) 

where Vj is the potential difference between states j + f 
and j while u)j{0, 0) corresponds to a spontaneous transi- 
tion probability from state j to state j + 1 in the absence 
of any potential difference and generalized force. Since 
the potential is periodic, the w^'s are subject to the con- 



dition V. 



0. In the absence of any potential dif- 



ference and generalized force, the particle can diffuse in 
both directions, so Uj(0,0) = Wj(0,Q) = uJj(0, 0). This 
choice corresponds to a FP equation in two dimensions 
where both the spatial and chemical coordinates, as well 
as their associated generalized forces, are treated on an 
equal footing^. 



According to eq. (|23|), all sums a+ + aj (m+ + mj) 
should be interpreted as the effective size of the spatial 
(chemical) substep taken by the motor protein. Some 
recent experiments [|l9| seem to suggest the existence of 
these small substeps. A substep can be of different types: 



only positional (m 



0), only chemical (a 



0) or mixed. The / and Afj, dependence in the transition 
rates, Uj's and Wj's, accounts also for more complicated 
schemes which cannot be ruled out, in principle. 

The substep size is unknown and it is related to the 
conformational changes involved after binding ATP, but 
on a pure theoretical basis, space may be discrctized so as 
to obtain the largest unit of substep such that all "natu- 
ral" substeps are multiple of this elementary unit. After 
all N substeps, a full spatial and chemical period has 
been covered, so that: 



**This equation may be written in the following form: 
dP(x,£,t) _ dJi(x,£,t) dJ 2 (x,£,t) 



dt 
Ji=D 1 

■h =D 2 



T dP p dv(x,p | ■ 

dx dx 



(49) 
(50) 
(51) 



where x and £ are, respectively, the spatial and chemical 
coordinates. 



4 + a i 



3 n' 



(52) 



where p is the typical spatial step performed by the mo- 
tor. We assume that the ujj for this elementary unit of 
substep do not depend on / and Afi. This hypothesis 
is commonly assumed in discrete models for substeps of 
any size [^Us|,[iO|] , whereas in our model it is assumed only 
for these elementary substeps. For a natural substep a 
more complicated force dependence of transition rates is 
possible, at least in principle. 

At the same time, since in absence of any internal po- 
tential or external generalized force the probability to 
proceed in one direction or the other is simply given by 
the probability to diffuse by 1/N of a full period in both 
coordinates, we assume that ujj does not depend on j and 
omit the subscript in the transition rates ujj. 

It is possible to show that such a model, in general, 
is not compatible with very simple requests, i.e. that 
the velocity of the motor saturates for high ATP concen- 
trations, and that the velocity obeys a Michaelis law, as 
follows from the data of |q,pj]]. Indeed, using eq. (|5^) it 
is possible to show that the velocity is given in general 
by (sec §): 



Jt 



puj I e 



/3A M 



1 



S 



(53) 



where 



JV 



/3[v n +(l-m+)Afi] 



71=1 

AT-1 



i=i j=i 



(54) 



in the case where / = 0. If all m„ > 0, in the large A/i 
limit the velocity is given by: 



Ji 



PLO 



y e P(v n -niZAn) 



(55) 



which is exponentially large in Afi unless all m+ = 0. 
But this condition, together with the condition that m+ + 
m n = lv implies that the velocity can be written in the 
form (q = e^): 



Jt. = 



A(q-l) 



K N q + K N ^iq « 



Kiq « 



, q = e 



13 An 



(56) 



where the constants A''s do not depend on q and are 
model dependent. Eq. ( J56] ) is not a Michaelis law, i.e. of 
the form J t = K ^+ n — B, even in the large q limit. This 

1. Therefore one 



K M +q 

is true for any value of N except N 



such discrete model is compatible with a Michaelis law 
only if we concentrate all the chemical reaction in a single 
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substep. Thus let us assume now that the chemical reac- 
tion occurs in a single step and without loss of generality 
we suppose that this ATP driven step is the first one in 
the cycle. We notice that this in turn would imply that 
the smallest substep is the one which occurs during ATP 
hydrolysis, which is reasonable^. 

A graphical representation of this class of discrete mod- 
els with only one chemical transition is given in fig. [l], for 
the case where N = 3. 



i - 2 



n - 1 




p + a p + 2a 



FIG. 1. Graphical representation of a discrete model, with 
N = 3 states along the spatial direction Xi. The chemical 
direction X2 represents the number of ATP molecules. Con- 
formational states of the motor protein are represented by 
the numbers 1,2,3. The system is periodic both in the spa- 
tial direction x\ (periodicity p; a = p/N — p/3) and in the 
chemical direction (periodicity 1). Only transition rates Ui 
and wi involve a chemical reaction, all the others are purely 
positional. 



Let A[i be apportioned over the forward and backward 



transitions, so that 



u X — e, m 1 — 
calculations it is possible to show that in this case we 
obtain for velocity an expression of the type: 



e. With some 



A(q-l) 



Kq + K\q e + K 2 q x 



(57) 



which is still incompatible with a Michaelis law unless 
e = or e = 1 . Choosing one or the other leads essentially 
to the same Michaelis law, except some multiplicative 
parameters which in turn depend on the products oje^ Vj , 
and ought to be determined by experiments. 

In the same way, the force apportionment can be in- 
vestigated. The problem is that at variance of chemical 
potential, experimental data arc obtained only for small 
forces. Therefore it is not possible to use the same criteria 



TT A net advancement of the protein center of mass is com- 
monly thought to occur upon release of the reaction products, 
whereas the hydrolysis reaction does not seem to imply any 
net macroscopic rearrangement on its own ||2Cf| . 



since, to our knowledge, for high values of force the ve- 
locity may grow up even to ±oo, meaning that the motor 
detaches from the fiber. Nevertheless, by the same line 
of reasoning, if we assume that the velocity should reach 
a constant value for high values of force in the positive 
direction, we find the condition: 



Vj = l, 



N. 



(58) 



On the contrary, if we assume that a constant velocity 
is reached only when the external force opposes the nat- 
ural direction of movement (which is what one expects 
on the basis of the available data jllj ) , then we find the 
condition that: 



Vj = l,...,iV 



(59) 



All intermediate cases, the symmetric one included, 
lead, therefore, to velocities growing up to ±oo as the 
force tends to ±oo. Remarkably, this is what happens 
when using a FP equation with very high forces, since 
the probability distribution is not sensitive to the po- 
tential shape and becomes flat, i.e. velocity is a linear 
function of force. Of course in a FP description every 
apportionment is totally equivalent to any other, as long 
as we consider the small a limit. This is not true in 
the limit for very high forces, for which our derivation of 
the continuum limit is no longer a good approximation. 
Therefore it is not surprising that in discrete models a 
particular choice of apportionment has dramatic conse- 
quences on the asymptotic behaviors of the quantities 
of interest, in the same way a different apportionment 
of the chemical potential leads to asymptotic behaviors 
which are not compatible with the expected Michaelis 
law unless, as shown above, the chemical reaction is con- 
centrated in a single step. In the following section we 
will study in detail some examples of models with spatial 
periodicity 2 and 3 highlighting their main predictions. 



IV. MODEL PREDICTIONS AND EXAMPLE 
STUDIES 

In all the following, lengths are measured in units of fil- 
ament periods p, while potential difference, Vj, and chem- 
ical potentials, A/i, are measured in units of kT, forces, 
/, in units of kT/p. 

According to our previous discussion, since experimen- 
tal data suggest that the velocity should reach a constant 
value for high negative forces, we assume the force ap- 
portionment to be asymmetric and concentrated in the 
forward transitions. In the same way we assume the first 
step to be chemically driven, so as to obtain a Michaelis 
law. We will show that these two hypotheses, left alone, 
are sufficient to obtain a force dependent Michaelis con- 
stant for velocity and rate of ATP consumption, and also 
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some interesting predictions about effective step-size and 
randomness factor. 
In these models: 



u n = u>e 



(60) 
(61) 



with the condition ^ n v n = 0. The expression for veloc- 
ity, using cqs. ( |60| ) and (|6l|), is: 



Ji 



V 0J (e/3( A ^+/) - 1) 



JV 



s 

«„+(l-5 n i)A^+(l-^)/ 



(62) 



n=l 



1 + 



i=l 



(63) 



Interestingly, the correct law for velocity is of the type: 
Aq 



Ji 



Km + q 



(64) 



where Km is the Michaelis constant. The same result 
has also been obtained in |(| in the context of continuous 
models. Eq. (|6^) may be used to calculate the 'stall' 
force, i.e. the force for which the velocity is zero. This is 
simply given by: 



f stall — — A/i. 



(65) 



This picture is, of course, very simplified. The exper- 
imental data reposted in show that the stall load 
(— /) increases with increasing ATP concentrations, but 
probably the dependence on A/z is not as simple as in 
eq. (|6^). However data in this parameter region are sub- 
ject to large experimental errors, due to rapid detach- 
ment of beads from the microtubule under stall condi- 
tions, so a comparison with experiment, at present, can 
be only of qualitative nature. 

The rate of ATP consumption follows from eq. p^ ) 
and can be written as: 



u\R\ — W1R2 

J2n R n 



(66) 



Using eq. (|6^) and eq. (|6(]), it is possible to show that 
the rate of ATP consumption and the velocity, when mea- 
sured in ATP molecules hydrolyzed per second and peri- 
ods per second respectively, are exactly the same quan- 
tity, so that the effective step size (the number of periods 
taken per hydrolyzed ATP molecule) is 1, which is consis- 
tent with experiments on kincsin |21|. This is no longer 
true if we use a more complicated scheme with pure ther- 
mal processes in addition to the normal ATP consuming 
ones, leading to transitions between different states: in 
this case the effective step size can be smaller than one. 



Using eq. (|6^) it is possible to obtain an explicit ex- 
pression for the force dependent Michaelis constant; for 
simplicity we concentrate on models with spatial peri- 
odicity 2 (2-modcls) and 3 (3-models). For 2-modcls 
V\ = v,V2 = —v, whereas for 3-models we assume 
V\ = v,V2 = —v/2,V3 = —v/2. This means that the 
chemical potential difference is used in the first transi- 
tion to overcome the internal potential barrier v, while 
the other transitions are favored by the internal potential 
shape. The Michaelis constant of eq. (p4) is given by: 



K M (2) 
Km (3) 



e 2v + 2e 
2ei v 



o -i- 
3e S 



e 3 



+2v 



1 + 2e7+i 



(67) 
(68) 



for 2- and 3-models respectively. In both examples, the 
Michaelis constant grows exponentially with — / (see dis- 
cussion of section III), in accordance with the recent ex- 
perimental observation of an increase in the Michaelis 
constant with applied load . Interestingly both mod- 
els predict a constant value for Km at high positive values 
of the force. We remark that experimental data on the 
Michaelis constant for positive forces, to our knowledge, 
are still unavailable. 

Recently developed experimental techniques allowed to 
measure the randomness parameter [ p2[ , defined as the 
long time limit of the ratio between the variance of the 
protein position on the filament, < x 2 (t) > — < x(t) > 2 , 
and the product of its average position, < x(t) > and 
periodicity p: 



lim 



< x 2 {t) > - < x(t) > 2 
< x(t) > p 



(69) 



In facts, the time resolution of experiments corre- 
sponds to the same order of magnitude of one hydrolysis 
event (typically milliseconds), while, at present, confor- 
mational changes leading to ATP hydrolysis cannot be 
directly measured. Nonetheless they affect the statistical 
distribution of the protein movement and the random- 
ness parameter U The macroscopic diffusion coefficient 
D is defined so that: < x 2 (t) > - < x(t) > 2 = 2Dt , 
while < x(t) >= Vt, where V is the velocity. Therefore: 



2D 
pV' 



(70) 



The expression for D is rather complicated, but still 
calculable for jump processes [Q. The randomness pa- 
rameter presents some interesting properties pl] , p2| , p3[ : 
first, r = for a perfectly clocklike motor, whereas 



^See [[l5| for a very interesting discussion on the relevance of 
the randomness parameter and the effective diffusion constant 
for mechanochemical transducers. 
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r = 1 for a 'Poisson' motor with one biochemical transi- 
tion and exponentially distributed time intervals between 
events. In general r _1 provides a continuous measure 
of the number of rate-limiting transitions in the over- 
all mechanochemical cycle (see |23j for a simple proof in 
the case of a jump process with equal forward and back- 
ward transition rates). In our models the expressions for 
the randomness parameter can be easily calculated, but 
their expressions are rather cumbersome. Some impor- 
tant features can be directly inferred and compared with 
experiments. 

At small loads — / and small A/x the velocity decreases 
linearly with both — / and A/x, so that the randomness 
factor from cq. (|7(i| ) should approach oo, assuming that 
the macroscopic diffusion coefficient approaches a con- 
stant value, in general different from zero. This should 
be also verified under stall conditions, when the random- 
ness factor again approaches oo. 



r 




0.4 



0.2 I , , , , , I 

5 10 15 20 25 30 

A^(kT) 

FIG. 2. Randomness parameter in a 2-model versus A/x, 
with v = 10 kT. Randomness approaches the value 1/2 when 
A/i ~ 2v, indicating that both transitions are rate limiting. 
At high A/i there is only one rate limiting transition, so that 
randomness approaches 1. The force / is measured in units 
of kT/p. 

At low values of A/i, but still sufficient to force the 
protein out of the stall condition, the randomness fac- 
tor exhibits different behaviors, depending on force, as 
shown in figures || and ||. When Afj, < v, the rate lim- 
iting step is essentially the first one for both models, so 
the randomness factor approaches 1 when the stall con- 
dition is overcome. At intermediate values, v < A/i < 2v 
for 2-models and v < A/i < |u for 3-modcls, the rate 
limiting transition is still the first one, but its rate limit- 
ing power is decreasing with respect to the others, due to 
the increase in A/i, so that the randomness parameter is 
decreasing, until at the border of this parameter region, 
A/i ~ v for 2- and A/i ~ Id for 3-models, all forward 
transition rates have essentially the same value: all steps 
are equally rate limiting, so that there are 2 rate limiting 



steps for 2-modcls and 3 for 3-models. This is evident 
also in the figures; in this part of the graphs, the ran- 
domness parameter approaches 1/2 and 1/3 for 2- and 
3-models respectively. 

At very high ATP concentrations, the first step has a 
high occurrence probability, whereas the other steps are 
rate limiting. This implies that the randomness parame- 
ter should approach a constant value 1 for 2-models and 
1/2 for 3-models. All these behaviors may be observed 
in the figures, at least for small values of force. Under 
very high loads, the velocity of the motor protein is very 
low in a wide range of A/i values. This implies that the 
randomness factor is very high, until A/i is high enough 
to force the system out of the stall condition. A direct 
comparison with experiments would be desirable. The 
central part of our figures reproduces data obtained for 
randomness in Our predictions in this parameter 

region also agree with theoretical derivations for contin- 
uous two-state ratchet models on kinesin |23fl , but in the 
case of continuous models it is very difficult to calculate 
the randomness parameter at very high ATP concentra- 
tions, due to numerical problems. This is no longer true 
for discrete models, for which all complications are alge- 
braic, rather than numerical. The experimental measure 
of randomness under very small loads or for small positive 
values of force should also be useful to infer, on a quan- 
titative basis, the spatial periodicity necessary to fully 
characterize the system. Data in JllJ seem to suggest a 
3-model, since the randomness factor approaches a min- 
imum value close to 1/3, but measurements at smaller 
loads would be useful to confirm or reject this hypothe- 
sis. 





f=0 — 




f-5 




f -8 - 




f=-1 




f=-12- 









5 10 15~~ 20 25 



A/z(fcT) 

FIG. 3. Randomness parameter in a 3-model versus A/i, 
with v = 10 kT. Randomness approaches the value 1/3 when 
A/i ~ 3u/2, indicating that all transition rates are rate lim- 
iting. At high A/i, there are essentially two rate limiting 
transition rates, so that randomness approaches 1/2, at least 
at low loads. The force / is measured in units of kT/p. 
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V. CONCLUSIONS 

In this paper, we have introduced a master equation 
approach to describe interesting properties of molecular 
motors. Our approach is similar in spirit to a Kramcrs- 
Moyal expansion of the FP equation 0] , but it is specif- 
ically studied for the problem of motor proteins, where 
both the mechanical and chemical coordinates are im- 
portant. In this framework we have studied the general 
conditions needed to obtain a discrete chemical kinetics 
model from a continuous one, when all generalized coor- 
dinates are treated on an equal footing. In this limit, a 
FP equation is recovered, when the requirement of a gen- 
eralized detailed balance condition, cq. (|23|), is fulfilled in 
the out of equilibrium regime. This condition has been 
shown to imply the validity of the Onsager reciprocity re- 
lations for the periodic stationary solutions in the close to 
equilibrium regime. This property is not surprising, since 
both Onsager relations and detailed balance are supposed 
to descend from microscopic reversibility. Motor pro- 
teins operate far from equilibrium, and far from the linear 
regime, where reciprocity relations are expected to hold. 
Nevertheless reciprocity relations constitute a minimal 
requirement in a model for mechanochcmical transduc- 
tion, as firstly stated by Hill |25|] and, more recently, in 
the context of continuous models |3). Moreover, to our 
knowledge, a proof for the case of discrete models, has 
never been given. 

The continuum limit has been shown to be insuffi- 
cient to fix all parameters in a discrete chemical kinetics 
model for any cyclic out of equilibrium thermodynamic 
system. However, in the case of motor proteins, semi- 
phenomenological considerations led us to the formula- 
tion of discrete models which are compatible with the 
stochastic continuous ones in the continuum limit, sat- 
isfy the Onsager reciprocity relations and allow an easy 
comparison with experimental data. 

We believe this work to be useful to study the general 
conditions which ought to be verified by a discrete ki- 
netics model and, most importantly, to study the force 
dependence of transition rates in mechano-chcmical pro- 
cesses. This is of relevance not only for research on motor 
proteins, but also on other important biomolcculcs. 
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Onsager relations eq. ([28|) hold for any discrete model. 
We also provide a derivation of the Onsager coefficients 
and give their values for the case of 2- and 3-models^] 
studied in section IV . In the following we use the nota- 
tion d a to denote a partial derivative with respect to a 
generalized force f a . 

Differentiating both sides of eq. (|2^), we obtain the 
condition: 



(d a L XY (F))R Y (0) - (d a L YX (F))R x (0) 
L X Y(0)R Y (0)(x a - y a ). 



F=0 



(Al) 



Differentiat ing the stationarity condition, eq. (g2|) , and 
applying cq. (Al) we obtain, after some manipulation: 



[Y /Y L XY (F)d a Ry(F)] F=o = 

a •^a ) • 

If AxY = Ax+LN.Y+LN^L £ Z d \ 

J2 XY ^ AxY = J2 XY *r AxY = 



(A2) 



E 



Xx + Xy 



XY 



A 



(A3) 



as a consequence of eq. (fL2|). 

Replacing q x with x a in cq. ( |l9| ) and differentiating 
with respect to fp, we obtain: 



A 



a/3 



{x a - y a ) {dpL XY ) R Y + 



xy 4 

(dpL YX ) R x + L XY (d j3 R Y ) - L YX (dpR x ) 



f=o 
(A4) 



Using condition (Al), we are left with: 

X a f3 — 



Xy 



E 



XY 4 

Xx "h X\ 
XY 4 



(x a - y a ) {x/3 - yp) Lxy(0)R y (0) + 



{Xql ya) 



'XY 



{dpR Y y 



^YX 



(dpRx) 



F=0 



(A5) 



Applying the property (A3) to the quantities (x a — 

y a )LxY{0)(d(3RY) and (x a - y a )L YX (0)(dpR x ) the On- 
sager coefficients are: 



APPENDIX A: ONSAGER COEFFICIENTS AND 
RELATIONS 

In this appendix we will prove that, provided a general- 
ized detailed balance condition in the form (E3f) holds, the 



"In these simplified models J c — J t , but this does not imply 
that the Onsager reciprocity relations are trivially satisfied. 
If they are, J c = Jt simply implies that all coefficients should 
be equal. 
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A a/3 

= E 



Xx ' Xy 



XY 



(x a - y a ) {xp - yp) L X y(0)Ry{0) + 



1 



E vv ( x <* ~ X Y L X Y{d(3R Y ) 



2 ^XY 

Xx L Yx{dpRx 



(A6) 



A subsequent application of eq. ( A2) and eq. (J23|) leads 



to: 



A a/3 



Xx ' Xy 



XY 



(x a - y a ) (xp - yp) L X y{0)Ry(0) + 



2 ^XY 

Lxy 



X, 



Rx 



X Y {d a Rx){dpR,Y)~ 
Ry 

(d a R Y )(dpR x ) 

- F=0 



(A7) 



From eq. ( |2 3[ ) and another application of eq. (A3), we 
finally obtain: 



Aq/3 — 



Xx +Xi 



^y(Q) 
i?x(0) ' 



(0) 



(i« - y a ) {xp - yp) L X y{Q)Ry{0)+ 

(A8) 



with Axy(F) = d a R x dpR Y +d a R Y dpR x . Eq. (|A§) 
is evidently symmetric under a change a <-» /3, so that 
finally the Onsager relations, eq. (Eq), are verified. The 



quantities j4xf(0) can be obtained by solving eq. (13) at 
stationarity for a finite number of states. 



a. 2— models 



For the 2-models defined in section [V, the coefficient 



A12 can be easily calculated; it is obtained from eq. (AS), 
by making explicit the two contributions: 



A (1) - 

A 12 — 



E 



Xx ' X\ 



XY 



(x t - yt) (x c - y c ) L XY (0)Ry{0) 



XY 



Lxy(0) 
4 Rx(0) 



A XY (0), 



(A9) 
(A10) 



where the subscript t (c) means translational (chemical). 
After some manipulations, and bearing in mind that xt — 
yt represents the spatial displacement from state Y to 
state X, while x c — y c is the ATP consumption from 
state Y to state X (positive when ATP is consumed) , we 
obtain: 



x (i) _ twi 

A 12 — 



Ui + W 2 



X (2) _ 
A 12 — 



2 U\ + 1t 2 + Wi + w 2 

U 1 (u 2 W2 — UxWi) 
2(m + W2)(U! + U 2 + Wi + U> 2 ) ' 



where all it, and Vi are calculated at (/, A/i) 
finally: 



^12 — 



W^Wi + 2uiU>iW2 + U1U2W2 

2(ui + w 2 )(ui + it 2 + Wi + W2) 



Using cqs. (60) and (pl|), we find: 
\ _ ^ 

A\2 — —■ 



(AH) 
(A12) 
(0, 0) and 

(A13) 
(A14) 



The same Onsager coefficient may be obtained by lin- 
earizing the expression for velocity, cq. (p2n . 



b. 3— models 

The calculation for the Onsager coefficient in 3-models 
by a direct application of eq. ( AS ) is rather long. Its value 
may be obtained by a direct linearization of eq. (62): 



A12 — 



- v + 2e-»/ 2 + 2e»/ 2 + e v + 3 ' 



(A15) 
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